
!==================================================================
!Interface to cusolverDn and CUDA C functions
!==================================================================
    
  module cuda_cusolver

  interface

     ! cudaMalloc
     integer (c_int) function cudaMalloc ( buffer, size ) bind (C, name="cudaMalloc" ) 
       use iso_c_binding
       implicit none
       type (c_ptr)  :: buffer
       integer (c_size_t), value :: size
     end function cudaMalloc

     ! cudaMemcpy 
     integer (c_int) function cudaMemcpy ( dst, src, count, kind ) bind (C, name="cudaMemcpy" )
       ! note: cudaMemcpyHostToDevice = 1
       ! note: cudaMemcpyDeviceToHost = 2
       use iso_c_binding
       type (C_PTR), value :: dst, src
       integer (c_size_t), value :: count, kind
     end function cudaMemcpy

     ! cudaFree
     integer (c_int) function cudaFree(buffer)  bind(C, name="cudaFree")
       use iso_c_binding
       implicit none
       type (C_PTR), value :: buffer
     end function cudaFree

     integer (c_int) function cudaMemGetInfo(fre, tot) bind(C, name="cudaMemGetInfo")
       use iso_c_binding
       implicit none
       type(c_ptr),value :: fre
       type(c_ptr),value :: tot
     end function cudaMemGetInfo


     integer(c_int) function cusolverDnCreate(cusolver_Hndl) bind(C,name="cusolverDnCreate")
     
     use iso_c_binding
     implicit none
     
     type(c_ptr)::cusolver_Hndl
     
     end function
     
     integer(c_int) function cusolverDnDestroy(cusolver_Hndl) bind(C,name="cusolverDnDestroy")
     
     use iso_c_binding
     implicit none
     
     type(c_ptr),value::cusolver_Hndl
     
     end function

    integer(c_int) function cusolverDnSgetrf_bufferSize(cusolver_Hndl,m,n,d_A,lda,Lwork) bind(C,name="cusolverDnSgetrf_bufferSize") 
      
      use iso_c_binding
      implicit none
    
      type(c_ptr),value::cusolver_Hndl
      integer(c_int),value::m
      integer(c_int),value::n
      type(c_ptr),value::d_A
      integer(c_int),value::lda
      type(c_ptr),value::Lwork
    end function

     integer(c_int) function cusolverDnDgetrf(cusolver_Hndl,m,n,d_A,lda,d_WS,d_Ipiv,d_devInfo) bind(C, name="cusolverDnDgetrf")

       use iso_c_binding
       implicit none
       type(c_ptr),value::cusolver_Hndl
       integer(c_int),value::m
       integer(c_int),value::n
       type(c_ptr),value::d_A
       integer(c_int),value::lda
       type(c_ptr),value::d_WS
       type(c_ptr),value::d_Ipiv
       type(c_ptr),value::d_devInfo
       
     end function 

     integer (c_int) function cusolverDnDgetrs(cusolver_Hndl,trans,n,nrhs,d_A,lda,d_Ipiv,d_B,ldb,d_devInfo) bind(C, name="cusolverDnDgetrs")

       use iso_c_binding
       implicit none
       type(c_ptr),value::cusolver_Hndl
       integer(c_int), value::trans
       integer(c_int), value::n
       integer(c_int), value::nrhs
       type(c_ptr),value::d_A
       integer(c_int), value::lda    
       type(c_ptr),value::d_Ipiv
       type(c_ptr),value::d_B
       integer(c_int),value::ldb
       type(c_ptr),value::d_devInfo
       
     end function
   
  end interface 
  
end module 

module cuda_getrs
    implicit none
	
	contains
  subroutine getrsGPU(A,B,ipiv,matrixSize,dim_change)
! The purpose of this routine is to provide for GPU based solution of 
! large dense systems of equations in legacy FORTRAN Applications.

! This is the development version of the routine which does not yet work ...

!cudaGetDeviceCount(nDevices)
  
  use iso_c_binding
  use cuda_cusolver
  
  ! ------------- matrix definition & host CPU storage variables 
  
  integer(c_int) m,n,lda,ldb
  
  real(8):: A(matrixSize,matrixSize)!,ATest(:,:)
  real(8):: B(matrixSize)
  integer::Ipiv(matrixSize),matrixSize,i,dim_change


!-------------------- CPU equivalents of device variables
  integer devInfo
  integer(c_int) nrhs
   integer*8 A_size,B_size,ipiv_size
  
! -------------------- handle to device
  
  type(c_ptr) :: cusolver_Hndl

  
! -------------------- pointers to device memory
  
  type(c_ptr),save::d_A
  type(c_ptr)::d_B
  type(c_ptr)::d_Ipiv
  type(c_ptr)::d_devInfo


! --------------------- function result variables
  
  integer cusolver_stat
  integer A_mem_stat
  integer B_mem_stat
  integer Ipiv_mem_stat

  
! -------------------- pointers to host CPU memory  
  
  type(c_ptr)::CPU_A_ptr,CPU_B_ptr
!  type(c_ptr)::CPUATest_ptr
  type(c_ptr)::CPU_ipiv_ptr
    
  target A,B,Ipiv !,ATest
  
  type(c_ptr)::cpfre,cptot
  integer*8,target::free,total
  integer res
  integer*8 cudaMemcpyDeviceToHost, cudaMemcpyHostToDevice
  integer*4 CUBLAS_OP_N, CUBLAS_OP_T
  parameter (cudaMemcpyHostToDevice=1)
  parameter (cudaMemcpyDeviceToHost=2)
  parameter (CUBLAS_OP_N=0)
  parameter (CUBLAS_OP_T=1)
  
  ! ================================================
  free = 0
  total = 0
  res = 1
  cpfre = c_loc(free)
  cptot = c_loc(total)
  
  res = cudaMemGetInfo(cpfre,cptot)
  if (res .ne. 0 ) then
      write (*,*)
      write (*, '(A, I2)') " cudaMemGetInfo error: ", res
      write (*,*)
      stop
  end if
  write (*, '(A, I12)') "  free mem: ", free
  write (*, '(A, I12)') " total mem: ", total
  
       
  lda = matrixSize
  ldb = matrixSize
  m=matrixSize
  n=matrixSize
  nrhs = 1
    
 
!  allocate(Lwork)

  

  
  !  "QR FACTORIZATION DENSE LINEAR SOLVER" 
  
  ! -------------- DEFINE [A] AND [B] --------------------------
  
  !           [A]         [x]  =    [B]
  !    | 1.0  2.0  3.0 | |1.0|    | 6.0|
  !    | 4.0  5.0  6.0 | |1.0| =  |15.0|
  !    | 2.0  1.0  1.0 | |1.0|    | 4.0|
  
 
      
  CPU_A_ptr = C_LOC(A)
  CPU_B_ptr = C_LOC(B)
   CPU_ipiv_ptr = C_LOC(ipiv)
 
  A_size = SIZEOF(A)
  B_size = sizeof(B)
  Ipiv_size = SIZEOF(Ipiv)
  !PRINT *,A_size
  
! Step 1: Create cudense handle ---------------
  cusolver_stat = cusolverDnCreate(cusolver_Hndl)  
  if (cusolver_stat .ne. 0 ) then
      write (*,*)
      write (*, '(A, I2)') " cusolverDnCreate error: ", cusolver_stat
      write (*,*)
      stop
  end if
  
! Step 2: copy A and B to Device
  if(dim_change==1) then
      A_mem_Stat    = cudafree(d_A)
      A_mem_stat    = cudaMalloc(d_A,A_size)
      if (A_mem_stat .ne. 0 ) then
          write (*,*)
          write (*, '(A, I2)') " cudaMalloc 1 error: ", A_mem_stat
          write (*,*)
          stop
      end if  
  end if
  B_mem_stat    = cudaMalloc(d_B,B_size)  
  if (B_mem_stat .ne. 0 ) then
      write (*,*)
      write (*, '(A, I2)') " cudaMalloc 2 error: ", B_mem_stat
      write (*,*)
      stop
  end if
  
! ---------- also allocate space for other device based variables  
  
  Ipiv_mem_stat = cudaMalloc(d_Ipiv,Ipiv_size)
  if (Ipiv_mem_stat .ne. 0 ) then
      write (*,*)
      write (*, '(A, I2)') " cudaMalloc 3 error: ", Ipiv_mem_stat
      write (*,*)
      stop
  end if
    
  ! ---------- copy A and B to Device
  if(dim_change == 1) then
  A_mem_stat = cudaMemcpy(d_A,CPU_A_ptr,A_size,cudaMemcpyHostToDevice)
      if (A_mem_stat .ne. 0 ) then
          write (*,*)
          write (*, '(A, I2)') " cudaMemcpy 1 error: ", A_mem_stat
          write (*,*)
    !      stop
      end if
  end if 
  B_mem_stat = cudaMemcpy(d_B,CPU_B_ptr,B_size,cudaMemcpyHostToDevice)
  if (B_mem_stat .ne. 0 ) then
      write (*,*)
      write (*, '(A, I2)') " cudaMemcpy 2 error: ", B_mem_stat
      write (*,*)
!      stop
  end if
   ipiv_mem_stat = cudaMemcpy(d_ipiv,CPU_ipiv_ptr,ipiv_size,cudaMemcpyHostToDevice)
  if (ipiv_mem_stat .ne. 0 ) then
      write (*,*)
      write (*, '(A, I2)') " cudaMemcpy 3 error: ", ipiv_mem_stat
      write (*,*)
!      stop
  end if
! Step 3: query working space of Sgetrf (and allocate memory on device)
  
 ! Lwork = 5
 ! cusolver_stat =  cusolverDnSgetrf_bufferSize(cusolver_Hndl,m,n,d_A,lda,CPU_Lwork_ptr) 
 ! if (cusolver_stat .ne. 0 ) then
 !     write (*,*)
 !     write (*, '(A, I2)') " DnSgetrf_bufferSize error: ", cusolver_stat
 !     write (*,*)
!      stop
!  end if
 

 
 ! Workspace = 4*Lwork
 ! WS_mem_stat = cudaMalloc(d_WS,Workspace)
 ! if (WS_mem_stat .ne. 0 ) then
 !     write (*,*)
 !     write (*, '(A, I2)') " cudaMalloc 6 error: ", WS_mem_stat
 !     write (*,*)
!      stop
 ! end if
  
! Step 4: compute LU factorization of [A] 
    

! Step 5: compute solution vector [X] for Right hand side [B]
  
  cusolver_stat = cusolverDnDgetrs(cusolver_Hndl,CUBLAS_OP_N,n,nrhs,d_A,lda,d_Ipiv,d_B,ldb,d_devInfo)  
  if (cusolver_stat .ne. 0 ) then
      write (*,*)
      write (*, '(A, I2)') " cusolverDnSgetrs error: ", cusolver_stat
      write (*,*)
!      stop
  end if
  
! Step 6: copy solution vector stored in [B] on device into [X] vector on host
  !PRINT *,B_SIZE
  B_mem_stat = cudaMemcpy(CPU_B_ptr,d_B,B_size,cudaMemcpyDeviceToHost)  
  if (B_mem_stat .ne. 0 ) then
      write (*,*)
      write (*, '(A, I2)') " cudaMemcpy 4 error: ", B_mem_stat
      write (*,*)
!      stop
  end if
  
  !do i = 1, n
   !  print *, B(i)
 ! enddo
  
! step 7: free memory on device and release CPU-side resources
  
  !A_mem_Stat    = cudafree(d_A)
  B_mem_Stat    = cudafree(d_B)
  Ipiv_mem_stat = cudafree(d_Ipiv)
  
  cusolver_stat = cusolverDnDestroy(cusolver_Hndl)
  
  ! Step 8: deallocate memory on host before exit
  
!  deallocate(A)


  end subroutine 
end module